function [rs]=trigonometric_polynomial_roots(S,C)
% solve a trigonometric polynomial of the following shape:
% \sum_{k=0}^{N} { sin(k-1)\theta S(k)+cos(k-1) \theta C(k)}=0
%   The trick is to turn the trigonometric polynomial into normal
%   polynomials, see the tech report for more details
%
%       lwei@my.fit.edu
% 
if(size(C,2)~=1||size(S,2)~=1)
    error('C and S must be column vectors');
end
if(numel(S)~=numel(C))
    error('S and C must have the same number of elements');
end
n=numel(S);
C=C.*1i;
T1=C-S;
T2=C+S;
T1=flipdim(T1,1);
T=zeros(2*n-3,1);
if(n>=4)
    T(1:n-3)=T1(1:n-3);
    T(n+1:end)=T2(n-1:end);
else
    error('n<4 is not supported yet');
end
T(n-2)=(S(1)-S(3)+C(3)+C(1));
T(n-1)=(-2*C(2));
T(n)=(S(3)-S(1)+C(1)+C(3));
[t]=roots(flipdim(T,1));

%% verify the solution
% for ii=1:numel(t)
%     ttt=t(ii);
%     ss=0;
%     for kk=1:n
%         ss=ss+S(kk)*(ttt^(kk-2)-ttt^(-(kk-2)))+C(kk)*(ttt^(kk-2)+ttt^(-(kk-2)));
%     end
%     disp(ss);
% end
%%
magt=t.*conj(t);
I=abs(magt-1)<1e-4;%throw the roots that does not satisify the trigonometry relationship
t=t(I);

if(numel(t)<1)
    rs=[];
    return;
end

%%
sinalpha=real((t-1./t)./2i);
cosalpha=real((t+1./t)./2);
rs=atan2(sinalpha,cosalpha);

%%

end